// Outcome model N(mu, sigma^2) (use SD rather than Var)
for (i in 1:Ni) {
Y_ijk[i] ~ bernoulli_logit(mu[i]);
}
}"
LogitModelle <- stan(model_code=modelle,init=0,data=daten,algorithm ="NUTS",iter=1000,warmup=750,chains=1)
test$lev3ForLev2 <- ave( test$hs_land, test$hs_land, FUN=seq_along)
daten <- list(Ni=2524,Nj=174,Nk=16,level2=test$hs_id,level=test$hs_land,
Y_ijk=test$small,X_1ijk=test$Sozial,lev3ForLev2=lev3ForLev2)
daten <- list(Ni=2524,Nj=174,Nk=16,level2=test$hs_id,level=test$hs_land,
Y_ijk=test$small,X_1ijk=test$Sozial,lev3ForLev2=test$lev3ForLev2)
modelle="data {
// Define variables in data
int<lower=0> Ni;       // Number of level-1 observations (an integer)
int<lower=0> Nj;       // Number of level-2 clusters
int<lower=0> Nk;       // Number of level-3 clusters
int<lower=1> level2[Ni];       // Cluster IDs
int<lower=1> level3[Ni];
int<lower=1> lev3ForLev2[Nj];
real Y_ijk[Ni];              // Logistic Outcome
real X_1ijk[Ni];             // Continuous predictor
}
parameters {
// Define parameters to estimate
real beta_0;                 // Population intercept (a real number)
real beta_1;                 // Population slope
real u_0jk[Nj];              // Level-2 random effect
real<lower=0> sigma_u0jk;
real u_0k[Nk];               // Level-3 random effect
real<lower=0> sigma_u0k;
}
transformed parameters {
real beta_0jk[Nj];         // Varying intercepts
real beta_0k[Nk];
real mu[Ni];              // Individual mean
// Varying intercepts definition
// Level-3 (10 level-3 random intercepts)
for (k in 1:Nk) {
beta_0k[k] <- beta_0 + u_0k[k];
}
// Level-2 (100 level-2 random intercepts)
for (j in 1:Nj) {
beta_0jk[j] <- beta_0k[lev3ForLev2[j]] + u_0jk[j];
}
// Individual mean
for (i in 1:Ni) {
mu[i] <- beta_0jk[level2[i]] + beta_1 * X_1ijk[i];
}
}
model {
// Prior part of Bayesian inference
// Flat prior for mu (no need to specify if non-informative)
// Random effects distribution
u_0k  ~ normal(0, sigma_u0k);
u_0jk ~ normal(0, sigma_u0jk);
// Likelihood part of Bayesian inference
// Outcome model N(mu, sigma^2) (use SD rather than Var)
for (i in 1:Ni) {
Y_ijk[i] ~ bernoulli_logit(mu[i]);
}
}"
LogitModelle <- stan(model_code=modelle,init=0,data=daten,algorithm ="NUTS",iter=1000,warmup=750,chains=1)
d <- transform(m, Cluster_ID = as.numeric(interaction(hs_land, hs_id, drop=TRUE)))
m <- cbind(test$hs_land,test$hs_id)
d <- transform(m, Cluster_ID = as.numeric(interaction(hs_land, hs_id, drop=TRUE)))
m
d <- transform(m, Cluster_ID = as.numeric(interaction(m[,1], m[,2], drop=TRUE)))
d
d[,3]
daten <- list(Ni=2524,Nj=174,Nk=16,level2=test$hs_id,level=test$hs_land,
Y_ijk=test$small,X_1ijk=test$Sozial,lev3ForLev2=d[,3])
modelle="data {
// Define variables in data
int<lower=0> Ni;       // Number of level-1 observations (an integer)
int<lower=0> Nj;       // Number of level-2 clusters
int<lower=0> Nk;       // Number of level-3 clusters
int<lower=1> level2[Ni];       // Cluster IDs
int<lower=1> level3[Ni];
int<lower=1> lev3ForLev2[Nj];
real Y_ijk[Ni];              // Logistic Outcome
real X_1ijk[Ni];             // Continuous predictor
}
parameters {
// Define parameters to estimate
real beta_0;                 // Population intercept (a real number)
real beta_1;                 // Population slope
real u_0jk[Nj];              // Level-2 random effect
real<lower=0> sigma_u0jk;
real u_0k[Nk];               // Level-3 random effect
real<lower=0> sigma_u0k;
}
transformed parameters {
real beta_0jk[Nj];         // Varying intercepts
real beta_0k[Nk];
real mu[Ni];              // Individual mean
// Varying intercepts definition
// Level-3 (10 level-3 random intercepts)
for (k in 1:Nk) {
beta_0k[k] <- beta_0 + u_0k[k];
}
// Level-2 (100 level-2 random intercepts)
for (j in 1:Nj) {
beta_0jk[j] <- beta_0k[lev3ForLev2[j]] + u_0jk[j];
}
// Individual mean
for (i in 1:Ni) {
mu[i] <- beta_0jk[level2[i]] + beta_1 * X_1ijk[i];
}
}
model {
// Prior part of Bayesian inference
// Flat prior for mu (no need to specify if non-informative)
// Random effects distribution
u_0k  ~ normal(0, sigma_u0k);
u_0jk ~ normal(0, sigma_u0jk);
// Likelihood part of Bayesian inference
// Outcome model N(mu, sigma^2) (use SD rather than Var)
for (i in 1:Ni) {
Y_ijk[i] ~ bernoulli_logit(mu[i]);
}
}"
LogitModelle <- stan(model_code=modelle,init=0,data=daten,algorithm ="NUTS",iter=1000,warmup=750,chains=1)
modelle="data {
// Define variables in data
int<lower=0> Ni;       // Number of level-1 observations (an integer)
int<lower=0> Nj;       // Number of level-2 clusters
int<lower=0> Nk;       // Number of level-3 clusters
int<lower=1> level2[Ni];       // Cluster IDs
int<lower=1> level3[Ni];
int<lower=1> lev3ForLev2[Nj];
real Y_ijk[Ni];              // Logistic Outcome
real X_1ijk[Ni];             // Continuous predictor
}
parameters {
// Define parameters to estimate
real beta_0;                 // Population intercept (a real number)
real beta_1;                 // Population slope
real u_0jk[Nj];              // Level-2 random effect
real<lower=0> sigma_u0jk;
real u_0k[Nk];               // Level-3 random effect
real<lower=0> sigma_u0k;
}
transformed parameters {
real beta_0jk[Nj];         // Varying intercepts
real beta_0k[Nk];
real mu[Ni];              // Individual mean
// Varying intercepts definition
// Level-3 (10 level-3 random intercepts)
for (k in 1:Nk) {
beta_0k[k] = beta_0 + u_0k[k];
}
// Level-2 (100 level-2 random intercepts)
for (j in 1:Nj) {
beta_0jk[j] = beta_0k[lev3ForLev2[j]] + u_0jk[j];
}
// Individual mean
for (i in 1:Ni) {
mu[i] = beta_0jk[level2[i]] + beta_1 * X_1ijk[i];
}
}
model {
// Prior part of Bayesian inference
// Flat prior for mu (no need to specify if non-informative)
// Random effects distribution
u_0k  ~ normal(0, sigma_u0k);
u_0jk ~ normal(0, sigma_u0jk);
// Likelihood part of Bayesian inference
// Outcome model N(mu, sigma^2) (use SD rather than Var)
for (i in 1:Ni) {
Y_ijk[i] ~ bernoulli_logit(mu[i]);
}
}"
LogitModelle <- stan(model_code=modelle,init=0,data=daten,algorithm ="NUTS",iter=1000,warmup=750,chains=1)
table(d[,3])
d[,3]
daten <- list(Ni=2524,Nj=174,Nk=16,level2=test$hs_id,level=test$hs_land,
Y_ijk=test$small,X_1ijk=test$Sozial,lev3ForLev2=d[,3])
modelle="data {
// Define variables in data
int<lower=0> Ni;       // Number of level-1 observations (an integer)
int<lower=0> Nj;       // Number of level-2 clusters
//     int<lower=0> Nk;       // Number of level-3 clusters
int<lower=1> level2[Ni];       // Cluster IDs
//     int<lower=1> level3[Ni];
//     int<lower=1> lev3ForLev2[Nj];
real<lower=0,upper=1> Y_ijk[Ni];     // Logistic Outcome
real<lower=0,upper=1> X_1ijk[Ni];    // Continuous predictor
}
parameters {
real beta_0;                 // Population intercept (a real number)
real beta_1;                 // Population slope
real u_0jk[Nj];              // Level-2 random effect
real<lower=0> sigma_u0jk;
//     real u_0k[Nk];               // Level-3 random effect
//     real<lower=0> sigma_u0k;
}
transformed parameters {
real beta_0jk[Nj];   // Varying intercepts
//     real beta_0k[Nk];
real mu[Ni];         // Individual mean
// Varying intercepts definition
// Level-3 (10 level-3 random intercepts)
//     for (k in 1:Nk) {
//       beta_0k[k] = beta_0 + u_0k[k];
//     }
// Level-2 (100 level-2 random intercepts)
//     for (j in 1:Nj) {
//       beta_0jk[j] = beta_0k[lev3ForLev2[j]] + u_0jk[j];
//     }
// Individual mean
for (i in 1:Ni) {
mu[i] = beta_0jk[level2[i]] + beta_1 * X_1ijk[i];
}
}
model {
// Random effects distribution
//     u_0k  ~ normal(0, sigma_u0k);
u_0jk ~ normal(0, sigma_u0jk);
for (i in 1:Ni) {
Y_ijk[i] ~ bernoulli_logit(mu[i]);
}
}"
LogitModelle <- stan(model_code=modelle,init=0,data=daten,algorithm ="NUTS",iter=1000,warmup=750,chains=1)
install.packages("coin")
library(fda)
annualprec = log10(apply(daily$precav,2,sum))
tempbasis =create.fourier.basis(c(0,365),65)
tempSmooth=smooth.basis(day.5,daily$tempav,tempbasis)
tempfd    =tempSmooth$fd
templist      = vector("list",2)
templist[[1]] = rep(1,35)
templist[[2]] = tempfd
conbasis = create.constant.basis(c(0,365))
betabasis = create.fourier.basis(c(0,365),5)
betalist = vector("list",2)
betalist[[1]] = conbasis
betalist[[2]] = betabasis
fRegressList = fRegress(annualprec,templist,betalist)
install.packages("yaml")
0.005<=0.001
1/100
#--- 1. Preparation -----------------------------------------------------------------------#
install.packages(c("coda","rstan","MCMCvis","dplyr"))
library(coda); library(rstan); library(MCMCvis); library(dplyr)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
set.seed(123454321)
#--- 2. Bayesian Statistical Modelling via RStan ------------------------------------------#
THKS <- read.csv("UseR_Ruhr_Daten.csv")
THKS <- subset(THKS, select=c("THKSbin", "PreTHKS", "CC", "TV", "CC.TV"))
#--- 2.1. Bayesian Logistic Regression: An Inefficient Example ---------#
Y <- THKS[,1]     # Response Variable
X <- THKS[,2:5]   # Predictor Variables
N <- nrow(X)      # Number of Respondents
D <- ncol(X)      # Number of Predictors
Input.Dataset <- list(Y=Y, X1=X[,1], X2=X[,2], X3=X[,3], X4=X[,4], N=N, D=D)
#-----------------------------------------------------------------------#
Logit_Model_1="
data {
//--- Definition of Model Information
int N;          // Index for Respondents
vector[N] X1;   // Vector for PreTHKS
vector[N] X2;   // Vector for CC
vector[N] X3;   // Vector for TV
vector[N] X4;   // Vector for CC.TV
int Y[N];       // Response Variable
}
parameters {
//--- Population-Level Effects
real alpha;     // Intercept
vector[4] beta; // Coefficients
}
model {
//--- Prior Distribution of Bayesian Inference
alpha ~ normal(0,10);     // Weakly-Informative Prior for Intercept
beta[1] ~ normal(0,10);   // Weakly-Informative Prior for Coefficients
beta[2] ~ normal(0,10);
beta[3] ~ normal(0,10);
beta[4] ~ normal(0,10);
//--- Likelihood Function of Bayesian Inference
Y ~ bernoulli_logit(alpha + beta[1]*X1 + beta[2]*X2 + beta[3]*X3 + beta[4]*X4);
}"
#-----------------------------------------------------------------------#
LogitModel1 <- stan(model_code=Logit_Model_1, model_name="Logit Model 1",
data=Input.Dataset, algorithm ="NUTS", init="random",
iter=5000, warmup=4000, chains=8, pars=c("alpha", "beta"))
THKS <- read.csv("UseR_Ruhr_Daten.csv")
setwd("~/ownCloud/rstan_material")
THKS <- read.csv("UseR_Ruhr_Daten.csv")
#--- 2. Bayesian Statistical Modelling via RStan ------------------------------------------#
THKS <- read.csv("UseR_Ruhr_Daten.csv")
THKS <- subset(THKS, select=c("THKSbin", "PreTHKS", "CC", "TV", "CC.TV"))
#--- 2.1. Bayesian Logistic Regression: An Inefficient Example ---------#
Y <- THKS[,1]     # Response Variable
X <- THKS[,2:5]   # Predictor Variables
N <- nrow(X)      # Number of Respondents
D <- ncol(X)      # Number of Predictors
Input.Dataset <- list(Y=Y, X1=X[,1], X2=X[,2], X3=X[,3], X4=X[,4], N=N, D=D)
#-----------------------------------------------------------------------#
Logit_Model_1="
data {
//--- Definition of Model Information
int N;          // Index for Respondents
vector[N] X1;   // Vector for PreTHKS
vector[N] X2;   // Vector for CC
vector[N] X3;   // Vector for TV
vector[N] X4;   // Vector for CC.TV
int Y[N];       // Response Variable
}
parameters {
//--- Population-Level Effects
real alpha;     // Intercept
vector[4] beta; // Coefficients
}
model {
//--- Prior Distribution of Bayesian Inference
alpha ~ normal(0,10);     // Weakly-Informative Prior for Intercept
beta[1] ~ normal(0,10);   // Weakly-Informative Prior for Coefficients
beta[2] ~ normal(0,10);
beta[3] ~ normal(0,10);
beta[4] ~ normal(0,10);
//--- Likelihood Function of Bayesian Inference
Y ~ bernoulli_logit(alpha + beta[1]*X1 + beta[2]*X2 + beta[3]*X3 + beta[4]*X4);
}"
#-----------------------------------------------------------------------#
LogitModel1 <- stan(model_code=Logit_Model_1, model_name="Logit Model 1",
data=Input.Dataset, algorithm ="NUTS", init="random",
iter=5000, warmup=4000, chains=8, pars=c("alpha", "beta"))
#--- 2.3. Bayesian Logistic Regression: An Inefficient Example ---------#
Y <- THKS[,1]     # Response Variable
X <- THKS[,2:5]   # Predictor Variables
N <- nrow(X)      # Number of Respondents
D <- ncol(X)      # Number of Predictors
Input.Dataset <- list(Y=Y, X=X, N=N, D=D)
#-----------------------------------------------------------------------#
Logit_Model_2="
data {
int<lower=1,upper=1600> N;  // Optimized Index for Respondents
int<lower=1,upper=4> D;     // Bounded Number of Predictors
int<lower=0,upper=1> Y[N];  // Optimized Response Variable
matrix[N,D] X;              // Matrix of Predictors
}
parameters {
real alpha;       // Intercept
vector[D] beta;   // Coefficients
}
model {
// Prior Distributions
alpha ~ cauchy(0,10);       // Optimized Weakly-Informative Prior
for(d in 1:D){
beta[d] ~ cauchy(0,2.5);  // Optimized Weakly-Informative Priors
}
// Optimized Vectorized Log-Posterior Density
Y ~ bernoulli_logit(alpha + X*beta);
}"
#-----------------------------------------------------------------------#
# Frequentist Logistic Model -> Produces Initial Values for Coefficients
Initialize <- glm(THKSbin~PreTHKS+CC+TV+CC.TV, data=THKS, binomial(link="logit"))
Init.Value <- coefficients(Initialize)
LogitModel2 <- stan(model_code=Logit_Model_2, model_name="Logit Model 2",
data=Input.Dataset, algorithm ="NUTS", init=Init.Value,
iter=5000, warmup=4000, chains=8, pars=c("alpha", "beta"))
#--- 3.4. Computation Results for the 'Efficient Example' --------------#
# Chain 1: 8.11716 Seconds      Beta1:    4147 Eff. Size
# Chain 2: 7.88067 Seconds      Beta2:    3570 Eff. Size
# Chain 3: 8.18669 Seconds      Beta3:    3847 Eff. Size
# Chain 4: 8.07629 Seconds      Beta4:    3515 Eff. Size
# Chain 5: 8.07522 Seconds      Alpha:    3315 Eff. Size
# Chain 6: 8.05129 Seconds
# Chain 7: 7.88022 Seconds
# Chain 8: 8.4576 Seconds
#-------------------------------------------------------
# Average: 8.090642 Seconds      Average:  3679 Eff. Size
#
# => Optimized Priors + Optimized Notation + Likely Initial Values
# => Massive Decrease in Average Computation Time by 61.61 Percent
# => Small Decrease in Average Sampling Efficiency by 50 Samples
#    probably due to a overly conservative Cauchy-Prior for alpha.
# ====> This is efficient Bayesian Statistical Modelling!
THKS <- read.csv("UseR_Ruhr_Daten.csv")
THKS <- subset(THKS, select=c("THKSbin", "PreTHKS", "CC", "TV", "CC.TV","Class"))
THKS <- arrange(THKS,Class)             # Reorganizing Data Set for Optimized Computation
Design.Matrix <- model.matrix(object = ~ 1+PreTHKS+CC+TV+CC.TV, data=THKS) # Design Matrix
Input.Dataset <- with(THKS,list(Ni=nrow(THKS), Nj=length(unique(Class)), Lvl2_ID=Class,
Y=THKSbin, DesMat=Design.Matrix, p=ncol(Design.Matrix)))
#-----------------------------------------------------------------------#
Two_Level_Logit_Model="
data {
//--- Definition of Level Populations
int<lower=1> Ni;              // Number of Level-1 Individuals
int<lower=1> Nj;              // Number of Level-2 Clusters
//--- Definition of Design Structure
int<lower=1> p;               // Number of Population Parameters
real DesMat[Ni,p];            // Design Matrix
//--- Definition of Cluster-ID
int<lower=1> Lvl2_ID[Ni];
//--- Definition of Dependent Variables
int<lower=0,upper=1> Y[Ni];   // Response Variable
}
parameters {
//--- Population-Level Effects
real beta[p];                 // Coefficient Matrix
//--- Group-Level Effects
real u_0jk[Nj];               // Level-2 Cluster Effects
real<lower=0> sigma_u0jk;     // Level-2 Cluster Standard Deviation
//--- The lower variance bound allows for faster covergence.
}
transformed parameters {
//--- Definition of Group-Variable Intercepts
real beta_0jk[Nj];
//--- Individual Population Mean
real mu[Ni];
//--- Computation of Group-Variable Intercepts
for (j in 1:Nj) {
beta_0jk[j]  = beta[1] + u_0jk[j];
}
for (i in 1:Ni) {
mu[i] = beta_0jk[Lvl2_ID[i]] + DesMat[i,2]*beta[2] +
DesMat[i,3]*beta[3] + DesMat[i,4]*beta[4];
}
}
model {
//--- Prior Distribution of Bayesian Inference
for (i in 1:p){
beta[i] ~ cauchy(0, 2.5);         // Optimized Weakly-Informative Cauchy-Prior
}
u_0jk ~ normal(0, sigma_u0jk);
sigma_u0jk ~ student_t(3, 0, 10);   // Optimized Weakly-Informative Hierarchical Student Folded-Prior
//--- Likelihood Function of Bayesian Inference
for (i in 1:Ni) {
Y[i] ~ bernoulli_logit(mu[i]);
}
}"
#-----------------------------------------------------------------------#
LogitModel3 <- stan(model_code=Two_Level_Logit_Model, model_name="Two-Level Hierarchical Logit Model",
data=Input.Dataset, algorithm="NUTS", init=0, iter=10000, warmup=7500, chains=8,
thin=20, pars=c("beta","sigma_u0jk"))
#--- 3. Bayesian Statistical Multilevel Modelling via RStan -------------------------------#
THKS <- read.csv("UseR_Ruhr_Daten.csv")
THKS <- subset(THKS, select=c("THKSbin", "PreTHKS", "CC", "TV", "CC.TV","Class","School"))
#--- 3. Bayesian Statistical Multilevel Modelling via RStan -------------------------------#
THKS <- read.csv("UseR_Ruhr_Daten.csv")
THKS <- subset(THKS, select=c("THKSbin", "PreTHKS", "CC", "TV", "CC.TV","Class","School"))
#--- 3.2. Bayesian Logistic Multilevel Regression: Three Levels --------#
THKS <- arrange(THKS,School,Class) # Reorganizing Data Set for Optimized Computation
Cluster.Setup <- unique(THKS[c("Class","School")])[,"School"] # Generation of Nesting-ID for Level 2
Design.Matrix <- model.matrix(object = ~ 1+PreTHKS+CC+TV+CC.TV, data=THKS) # Design Matrix
Input.Dataset <- with(THKS,list(Ni=nrow(THKS), Nj=length(unique(Class)), Nk=length(unique(School)),
Lvl3_ID=School, Lvl2_ID=Class, Lvl3_Lvl2=Cluster.Setup, Y=THKSbin,
DesMat=Design.Matrix, p=ncol(Design.Matrix)))
#-----------------------------------------------------------------------#
Three_Level_Logit_Model="
data {
//--- Definition of Level Populations
int<lower=1> Ni;              // Number of Level-1 Individuals
int<lower=1> Nj;              // Number of Level-2 Clusters
int<lower=1> Nk;              // Number of Level-3 Clusters
//--- Definition of Design Structure
int<lower=1> p;               // Number of Population Parameters
real DesMat[Ni,p];            // Design Matrix
//--- Definition of Cluster-IDs
int<lower=1> Lvl3_ID[Ni];
int<lower=1> Lvl2_ID[Ni];
int<lower=1> Lvl3_Lvl2[Nj];
//--- Definition of Dependent Variables
int<lower=0,upper=1> Y[Ni];   // Response Variable
}
parameters {
//--- Population-Level Effects
real beta[p];                 // Coefficient Matrix
//--- Group-Level Effects
real u_0k[Nk];                // Level-3 Cluster Effects
real u_0jk[Nj];               // Level-2 Cluster Effects
real<lower=0> sigma_u0k;      // Level-3 Cluster Standard Deviation
real<lower=0> sigma_u0jk;     // Level-2 Cluster Standard Deviation
}
transformed parameters {
//--- Definition of Group-Variable Intercepts
real beta_0jk[Nj];
real beta_0k[Nk];
//--- Individual Population Mean
real mu[Ni];
//--- Computation of Group-Variable Intercepts
for (k in 1:Nk) {
beta_0k[k]  = beta[1] + u_0k[k];
}
for (j in 1:Nj) {
beta_0jk[j]  = beta_0k[Lvl3_Lvl2[j]] + u_0jk[j];
}
for (i in 1:Ni) {
mu[i] = beta_0jk[Lvl2_ID[i]] + DesMat[i,2]*beta[2] +
DesMat[i,3]*beta[3] + DesMat[i,4]*beta[4];
}
}
model {
//--- Prior Distribution of Bayesian Inference
for (i in 1:p){
beta[i] ~ cauchy(0, 2.5);
}
u_0k  ~ normal(0, sigma_u0k);
sigma_u0k ~ student_t(3, 0, 10);
u_0jk ~ normal(0, sigma_u0jk);
sigma_u0jk ~ student_t(3, 0, 10);
//--- Likelihood Function of Bayesian Inference
for (i in 1:Ni) {
Y[i] ~ bernoulli_logit(mu[i]);
}
}"
#-----------------------------------------------------------------------#
LogitModel4 <- stan(model_code=Three_Level_Logit_Model, model_name="Three-Level Hierarchical Logit Model",
data=Input.Dataset, algorithm="NUTS", init=0, iter=10000, warmup=7500, chains=8,
thin=20, pars=c("beta","sigma_u0jk","sigma_u0k"))
